"""
Purpose: Predicting drug-target interactions via the R package "BioMedR"

Author: H. Lin, PhD . Https://orcid.org/0000-0003-4060-7336

Version: Created on 14th Aug. 2020; updated on 21st Jan. 2024; Lasted updated on 21st Aug. 2025

Drug-Target predictions --- Computational predictive analyses on drug-target interactions
"""

# Load the library required. 
library(BioMedR) # Note that, multiple package-/version- dependency issues may exist according to different users' computing environment. E.g., the rlang and GOSemSim packs.

gpcr <- read.table(system.file('vignettedata/GRCR.csv', package= 'BioMedR') # Load the example dataset embedded in the package

protid <- unique(gpcr[,1]); drugid=unique(gpcr[,2]) # Protein and drug ID assignments to variable

protseq <- BMgetProtSeqKEGG(protid, parallel=5) # Retrieve protein seqeuence data from remote database, using IDs. Network connection is required

drugseq=c() # Create a NULL variable for storging SMILES data of drugs

for (id in 1:length(drugid)){drugseq[id]=BMgetDrugSmiKEGG(drugid[id])} # Remote retrieval of drug SMILES via drug IDs


"""  Alternatively, below two lines of codes could be used if your network connection is not good

protseq = readFASTA(system.file('vignettedata/GPCR_seq.fasta', package = 'BioMedR'))

drugseq = as.vector(read.table(system.file('vignettedata/GPCR_smi.txt', package = 'BioMedR'), col.names = 'SMILES'))
"""                


x0.prot = cbind(t(sapply(unlist(protseq), extrProtMoreauBroto)), t(sapply(unlist(protseq), extrProtCTDC))) # Combine two featuresets of protein sequences

x0.drug = cbind(extrDrugGraphComplete(readMolFromSmi(textConnection(drugseq))), extrDrugPubChemComplete(readMolFromSmi(textConnection(drugseq))))  
                
x.prot=matrix(NA, nrow = nrow(gpcr), ncol=ncol(x0.prot)) 

x.drug=matrix(NA, nrow = nrow(gpcr), ncol = ncol(x0.drug)) # Creating matrices via specifying their row number and column numbers

for(i in 1:nrow(gpcr)) x.prot[i,] = x0.prot[which(gpcr[,1][i] == protid),]

for(i in 1:nrow(gpcr)) x.drug[i,]=x0.drug(gpcr[,2][i] == drugid) # Assign values to matrices

y=as.factor(c(rep('1',nrow(gpcr)/2), rep('0', nrow(gpcr/2)))) # Create label set

x=getCPI(x.prot, x.drug, type='combine') # Generate drug-target interaction descriptors (combined descriptor matrix) using getCPI(). 

colnames(x) =paste('CCI', 1:dim(x)[2],sep='_') # Fill the column name with numbers

require(caret) # Load library. If not installed, should firstly install via command "install.packages("caret") ".  But note that, package-/version- dependent issues may exist. In my R environment, "rlang" and "vctrs" package version dependent issue happened. Hence extra procedures for installation or updating packs were required. 

x=x[,-nearZeroVar(x)]

# Machine Learning Training set splitting
set.seed(20180808) # Set.seed for result reproduction/regeneration

split_index=createDataPartition(y,p=0.75, list=F)  # createDataPartition ( ) function for splitting training and testing sets.  The parameter p's value indicate the percentage for training set. Here 75% were used for training set.

train_x=x[split_index,] # training set split for feature set
train_y=y[split_index] # training set split for label set
test_x=x[-split_index,] # testing set split for feature set
test_y=y[-split_index] # testing set split for label set


require(randomForest) # Load the library of randomforest

cv_result=rfcv(train_x,train_y, cv.fold=5, type='classification', tree=500, mtry=30) # Five-fold cross-validation with random forest classifier. Note that the correct function name should be 'rfcv( )' instead of rf.cv( ) written in the official tutorial document


# Train the random forest classifier
rf.fit=randomForest(x=train_x, y=train_y, ntree=500, mtry=30, importance=TRUE) 


# Predict on the test set (in fact the training set)
pre_res=predict(rf.fit, newdata=test_x, type='prob')[,2]

require(pROC) #  plot the Cross validation result and test result
require(RColorBrewer)
pal=brewer.pal(3, 'Set1')
opar<-par(no.readonly = TRUE)
par(mfrow=c(1,2))
plot.roc(train_y, cv_results$prob[,1], col=pal[2], grid=TRUE, print.auc=TRUE, main=' Cross Validation')
plotroc.(test_y, pre_res, col=pal[1], grid=T, print.auc=T, main='prediction')
par(opar)

                
       
            

""" References
> citation("caret") 

在出版物中使用程序包时引用‘caret’:

  Max Kuhn (2020). caret: Classification and
  Regression Training. R package version
  6.0-86.
  https://CRAN.R-project.org/package=caret

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {caret: Classification and Regression Training},
    author = {Max Kuhn},
    year = {2020},
    note = {R package version 6.0-86},
    url = {https://CRAN.R-project.org/package=caret},
  }

> citation("BioMedR")

在出版物中使用程序包时引用‘BioMedR’:

  Min-feng Zhu, Jie Dong and Dong-sheng Cao
  (2019). BioMedR: Generating Various
  Molecular Representations for Chemicals,
  Proteins, DNAs, RNAs and Their
  Interactions. R package version 1.2.1.
  https://CRAN.R-project.org/package=BioMedR

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {BioMedR: Generating Various Molecular Representations for Chemicals,
Proteins, DNAs, RNAs and Their Interactions},
    author = {Min-feng Zhu and Jie Dong and Dong-sheng Cao},
    year = {2019},
    note = {R package version 1.2.1},
    url = {https://CRAN.R-project.org/package=BioMedR},
  } 


Codes in this R script file was referred to and modified based on the official tutorial/manual of the BioMedR package: BioMedR: R/CRAN Package for generating various molecular representations for chemicals, proteins DNAs/RNAs and their interactions. The BioMedR manual was authored by Minfeng Zhu, Jie Dong, Dongsheng Cao, Package version: Release 1. 2019-07-03 



Codes were fully tested by H. Lin.

Parts of error and buggy codes from the original official manual of BioMedR were fixed and revised by H.Lin

The step-by-step code instructions/explanations/annotations in this script file were written by H. Lin.  

"""
